Onsager Loop-Transition and First Order Flux-Line Lattice Melting in High-T c 

Superconductors 



A. K. Nguyen and A. Sudb0 
Department of Physics 
Norwegian University of Science and Technology, N-7034 Trondheim, Norway 

Monte-Carlo simulations in conjunction with finite-size scaling analysis are used to investigate 
the (H , T)-phase diagram in uniaxial anisotropic high-T c superconductors, both in zero magnetic 
field (B = 0) and in intermediate magnetic fields (0 -C B <C B C 2) for various mass-anisotropies. The 
model we consider is the uniformly frustrated anisotropic Villain Model, which is dual to the Lattice 
London Model with an infinite London penetration length. The quantities we consider are various 
helicity moduli, the structure function, the specific heat, and the distribution of closed non-field 
induced vortex loops as a function of the loop-size. In zero magnetic field, and for all anisotropics 
considered, we find one single second order phase transition, mediated by an Onsager vortex-loop 
unbinding transition, or blowout. This is the superconductor-normal metal transition. A comparison 
with numerical simulations and a critical scaling analysis of the zero-field loop-transition yields the 
same exponent of the loop-distribution function at the critical point. In the intermediate magnetic 
field regime, we find two anomalies in the specific heat. The first anomaly at a temperature T m is 
associated with the melting transition of the flux-line lattice. The Lindemann-ratio at the melting 
is given by cl ~ 0.24. The second anomaly at a temperature T z is one where phase coherence 
in the BCS order parameter across the sample along the field direction is destroyed. We argue 
that Tm = T z in the thermodynamic and continuum limit. Hence, there is no regime where the 
flux-line lattice melts into a disentangled flux-line liquid. The loss of phase coherence parallel to 
the magnetic field in the sample is argued to be due to the proliferation of closed non-field induced 
vortex loops on the scale of the magnetic length in the problem, resulting in flux-line cutting and 
recombination. In the flux-line liquid phase, therefore, flux-lines appear no longer to be well defined 
entities. Above the melting temperature, the system always exhibits an incoherent vortex-liquid 
phase characterized by lack of phase coherence in the BCS order parameter parallel to the magnetic 
field. For increasing anisotropy, we resolve a delta-function peak in the specific heat. A finite-size 
scaling analysis of the delta-function peak specific heat anomaly at the melting transition is used 
to extract the discontinuity of the entropy at the melting transition. This entropy discontinuity is 
found to increase rapidly with mass-anisotropy, at least for not too layered compounds. 



I. INTRODUCTION 



A number of recent experiments have reported results 
of a first order melting transition of the Abrikosov flux- 
line lattice (FLL)El u. The reported magnitudes of the 
latent heat have all in general been in surprisingly good 
agreement with a prediction of Hetzel et al. for the dis- 
continuity in the entropy at the melting transition based 
on extensive Monte-Carla simulations of the uniformly 
frustrated 3D AF-modeM 

Schilling et al. have reported calorimetric measurements 
on an untwinned Y BaiCuzOi (YBCO) single crystal, in 
the intermediate field regime i?e[l — 7] Tesla, and find 
a FLL melting transition with a virtually field indepen- 
dent entropy jump AS ~ 0A5ks per vortex per layeru. 
In Bi 2 Sr2CaiCu20s (BSCCO) single crystal at very low 
magnetic inductions, Be[1— 375] Gauss (G), it was found 
that a FLL melting transition occurs with an enormous 
entropy jump AS(B=1G) ~ 6ks per vortex per layerQ. 
Furthermore, it was found that AS(Bl decreases for in- 
creasing B, and vanishes at B^375 Ga. It appears that 
AS(B) increases dramatically only when B — > 0, T — > T c . 
In fact, it might be argued on the basis of the data of Zcl- 
dov et al., that AS(B) diverges in this limit. 



Recently, Tesanovidlj and Nguyen et alH have proposed 
an explanation for the inordinately large entropy jump 
found in Ref.Q. The idea is that the FLL melting tran- 
sition takes place at roughly the same temperature as 
a "blowout" of non-field induced degrees of freedom in- 
volving closed vortex loops, resulting in a flux-line liquid 
phase with considerably larger entropy than what the 
field-induced vortices alone can provide. 
As a step towards understanding these experimental re- 
sults, we carry out extensive Monte Carlo simulations, to- 
gether with a finite-size scaling analysis, on the uniformly 
frustrated anisotropic Villain Model, to be defined below. 
This model will be argued to be appropriate for describ- 
ing the physics in extreme type-II superconductors such 
as the high-T c superconductors. Here, we present a short 
review of our results. 

For the case B = 0, we find that the Villain Model has one 
single second order phase transition of the 3DAV-type 
for all anisotropics considered. The phase transition in 
zero magnetic field is caused exclusively by a vortex loop 
"blowout" , to be explained below. This is confirmed by a 
detailed calculation of the distribution function for loops 
of a given perimeter, as a function of the perimeter, for 
various temperatures. In the low temperature regime, 
this distribution function is an exponentially decreasing 
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function of the perimeter, indicating that there exists a 
length scale in the problem associated with a typical size 
of thermally induced closed vortex-loops in the system. 
However, in zero magnetic field there exists a tempera- 
ture scale, which we denote as T c , at which the distri- 
bution function is an algebraically decreasing function of 
the perimeter of vortex loops, indicating that there no 
longer exists a length scale associated with typical sizes 
of thermally induced closed vortex-loops. In this case, 
such vortex-loops exists on all length scales in the prob- 
lem, up to and including the system size. 
This means that the system experiences a thermally in- 
duced proliferation of unbounded closed vortex-loops, a 
situation for mhich Onsager coined the term "vortex- 
loop blowout" E3. In zero magnetic field, such a blowout 
marks the transition from a parmal metal to a supercon- 
ducting state, or vice versaEjO. Many years ago, the 
vortex-loop blowout transition was suggested to .occur in 
the neutral super-fluid £/e 4 at the A-transitionE3. The 
loop-transition in He A has-more recently been reinves- 
tigated by several authorsEd. In the context of charged 
super-fluids in zero magnetic field a corresponding loop- 
transition was suggested to occur in isotropic lattice su- 
perconductor models several years ago by D.asgupta and 
HalperintS, and more recently by Korshunovli3. The sug- 
gestion that features of this zero-field transition may sur- 
vive in finite magnetic fields, and thus be of importance 
for the statistical mechanics of the vortex-liquid phase, 
has been suggested by TesanoviaI3, and considered re- 
cently by us in detailed, Monte-Carlo simulations of the 
Lattice London ModelB. 

The main purpose of the present paper is to study, via de- 
tailed Monte-Carlo simulations, the fate of the zero-field 
Onsager-Dasgupta-Halperin transition when a magnetic 
field is applied to an extreme type-II superconductor, us- 
ing the somewhat simpler more familiar uniformly frus- 
trated anisotropic Villain-Model, which is related to the 
Lattice London Model via a duality transformation. 
For finite fields B / 0, we find two sharp features in 
the specific heat and helicity modulus. In addition, 
we observe what appears to be a cross-over at a con- 
siderably \&l£££ temperature, in agreement with recent 
simulationstjO . 

The first sharp feature we find, at a temperature which 
we denote as T m , is identified as the first order melting 
transition of the FLL. The second sharp feature is more 
subtle. It takes place at a temperature which we denote 
as T z . At the temperature T z , we find that the phase- 
coherence in the BCS superconducting order-parameter 
across the sample in the direction of the magnetic field 
is destroyed. We also find, via our computations, that 
flux-line cutting and the amount of intersecting flux lines 
dramatically increases at T z . Consequently, above T z , 
phase coherence along the field direction is destroyed. 
Furthermore, for the anisotropics considered in this pa- 
per, it occures that T z —> T m from above as the system 
size is increased. T z never drops below T m , for reasons 
to be explained below. We emphasize that we at this 



stage are limiting this statement to the case of infinite 
penetration depth, since the results are obtained within 
the uniformly frustrated Villain-model only. 
From this we draw two conclusions. Firstly, when the 
flux-line lattice melts, it does not melt into a flux-line 
liquid which has phase-coherence along the direction of 
the applied magnetic field in any temperature regime. 
This however does not mean that the flux-line liquid is 
an entangled vortex system: At T z — > T m we find that 
flux-line cutting and intersectioning of flux lines with 
closed vortex-loops of diameter on the scale of the mag- 
netic length, increases abruptly. Hence, the flux lines 
cannot be considered as well defined in the liquid phase. 
Secondly, in the present simulations on the anisotropic, 
uniformly frustrated Villain-Model, we never observe an 
entangled flux-line lattice phase of the type described by 
Frey et alL3 for the layered case. However, we cannot 
access extreme anisotropies in our simulations, for rea- 
sons to be explained in Section II. D, and therefore do 
not rule out the existence of a "supersolid phase" such 
as proposed in Ref.EEl. 

We also observe a crossover that takes place at a temper- 
ature which marks the onset of strong diamagnetism, not 
associated with global phase-coherence in the supercon- 
ducting BCS order parameter, but with phase-coherence 
throughout finite domains. It takes place at a temper- 
ature well above both T m and T z , which we denote as 
Tbc2- Usually, this cross-over is identified with B C 2, the 
upper critical field, and signals the onset of strong dia- 
magnetic fluctuations. This crossover is the remnant of 
the zero-field second order phase-transition that marks 
the onset of the transition form metallic to supercon- 
ducting behavior. The filling fractions we consider in 
this paper, / = 1/32 and / = 1/72, may be converted 
into a magnetic field of the order of 10 Tesla, which is 
not particularly low. The present paper, therefore, does 
not address the issue of how the crossover at B C 2 in finite 
large fields evolves into the sharp second-order transition 
in zero field. This issue is of fundamental importance, 
and remains open. 

The rest of this paper is organized as follows. In Sec- 
tion II we describe the uniformly frustrated anisotropic 
Villain Model along with the approximations inherent in 
this description of a superconductor. We also give the 
connection between these models and the Lattice Lon- 
don Model. Then we define the physical quantities to be 
considered, and their measurements in the simulations. 
In Sections III and IV we show and discuss our results for 
the zero-field case and the finite-field case, respectively. 
Finally in Section V we summarize our main findings. 
The derivations of the helicity moduli we consider, both 
in terms of phase- variables and in terms of vorticities, are 
given in two appendices. 



II. THE MODEL AND DEFINITIONS 
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A. The Model 

Our starting point is the anisotropic lattice supercon- 
ductor model (LSM) (which semantically should be dis- 



tinguished from the Lattice London Model)0S@, 
fined by the partition function 
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Here, Jo is the energy scale for the system. Fur- 
thermore, a p is the anisotropy parameter along the /x- 
direction, and V denotes a lattice derivative. The vari- 
able 9{r) S [— 7r, 7r) is the phase of the complex supercon- 
ducting order parameter ^(r) at site f of a three dimen- 
sional numerical lattice with lattice constant d, m M (r) are 
integer variables defined on the directed link between site 
r and site r + e M , where is the primitive vectors for the 
cubic unit cell (|e M |=d, Fig. [l]). The contribution A^{f) 
to the gauge-invariant phase of the order-parameter is 
related to the vector potential A vp (f) by 
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where $ = 2.07- 10~ 15 Tm 2 is the flux quantum. Finally, 
A M is the London penetration depth along the /i-direction. 
In this model, we neglect fluctuations of the amplitude 
in the complex superconducting order parameter, i.e. 
$r(f)= \ty(f)\ e ie W -^ e ie(f?) . The Lattice London Model 
is obtainable from the Lattice Superconductor Model by 
explicitly performing the 9 and ^-integrations in Eq. 
([j]), as shown first by Korshunov and more recently by 
CarneiroBlii 

To study the physics of high-T c superconductors, we con- 
sider a three dimensional cubic lattice, with linear dimen- 
sion Ld, and with a uniaxial anisotropy T = A c /A a . In 
these simulations the applied magnetic field, and hence 
the net magnetic induction B, is taken along the crystal 
c-axis. Here, A a and A c are the penetration depths in 
the crystal ab-plane (CitC^-plane) and along the crystal 
c-axis, respectively. Subsequently, we will take the limit 
A a7 A c —> 00, but in such a way that the ratio A c /A a is 
maintained constant. We take our coordinate (x, y, z)- 
axis parallel to the crystal (a, 6, c)-axis, respectively. Pe- 
riodic boundary conditions in all direction are_assumed. 
The basic parameters of the LSM are given byEj, 
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Here, d may tentatively be interpreted as the distance be- 
tween two CuC^-layers in adjacent unit cells. The vortic- 
ities n x (r),n y (r) (corresponding to fundamental vortex- 
line segments parallel to the ab-plane, defined in Eq. ||) 



are assumed to exist in between CuC>2 double-or multiple- 
layers in compounds such as YBCO and BSCCO. We use 
the numerical lattice unit as a measure of the in-plane 
coherence length, £ a ~ d. Note that since the numerical 
factor relating £ a to d is not uniquely determined in our 
approach, the filling fraction / does not uniquely deter- 
mine the magnitude of the applied magnetic field. Vari- 
ation of / may thus be viewed as a variation of B, but 
alternatively a reduction of / may be viewed as an im- 
provement of the approximation to the continuum limit 
at fixed induction B.) 

As for the London Model, the LSM is appropriate for de- 
scribing the physics of extreme type-II superconductors 
(A a S> £a) in the field regime B <C B c2 , where B c2 is the 
upper critical magnetic field, implying a v ^> £ a . Thus, 
spatial variations of the amplitude of the superconducting 
order parameter may be neglected. In these simulations 
it is also postulated that the penetration length is es- 
sentially infinite, which in practice means that they are 
at least much larger than the average distance between 
vortex lines, when the field is finite. Hence we also have 
the requirement B 3> B c \. In terms of magnetic induc- 
tion, our simulations are thus strictly speaking limited 
to the field-regime B c \ <C B <C B c2 when finite fields 
are considered. For our zero-field results, the complete 
suppression of gauge-fluctuations implies that the pene- 
tration depth of the model must be at least larger than 
any system-sizes considered. 

The Monte Carlo simulation time tMC for the LSM on a 
cubic system with linear dimension L is of order L 6 . The 
suppression of the gauge-field fluctuations, using the limit 
A M =oo, reduces the required computer time dramatically, 
tMC to^L 3 . The neglect of gauge-fluctuations reduces 
the LSM— ta-the uniformly frustrated anisotropic Villain 
Mode£3ll3 : EJ, which is the model used in this paper. It 
is defined by the following partition function after per- 
forming the sum over m M (r) in Eq. |l| explicitly, 
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The advantage of the Villain-model model compared to 
the Lattice London Model used in earlier large-scale sim- 
ulations, we have performed on the Abrikosov vortex 
latticeEj, is that is allows considerably larger system sizes 
to be studied than with the Lattice London Model. The 
latter model has the intuitively appealing feature of al- 
lowing simulations on line-like objects, but as we have 
seen, the Villain-model and the Lattice London model 
are in principle equivalent representations of a lattice 
superconductor model. One other major advantage of 
the Villain-model compared to the Lattice London model 
is that, while it is straightforward to extract unambigu- 
ous information about vorticities from the phases of the 
Villain-model, it is impossible to reconstruct unambigu- 
ous phase-information from the vortex-degrees of free- 
dom in the Lattice London model. Thus the Villain- 
model straightforwardly provides information on vortic- 
ities as well as phase-coherence. Ultimately, the choice 
of model to use in simulations is dictated exclusively by 
convenience, and depends to a large degree on what prob- 
lems to consider. One problem where the Lattice Lon- 
don Model appears to present clear advantages over the 
Villain-model is the problem of numerical simulations of 
flux-creep in the presence of pinning. 
The uniformly frustrated anisotropic Villain Model is ap- 
propriate for describing the physics of extreme type-II su- 
perconductors in the limit where the penetration depth 
is larger or comparable to the system size (zero magnetic 
field) or when the penetration depth is much larger than 
the average distance between flux lines A S> a v (finite 
magnetic fields). 



Here, Ci is the closed path traced out by the links sur- 
rounding an arbitrary plaquette, and v represents the 
Cartesian components of the current in the directions of 
the links which comprise the closed path Ci. Further- 
more, n^(r) = 0, ±1 represents a vortex segment pene- 
trating the plaquette enclosed by the path Ci. The situa- 
tion is illustrated in Fig. Furthermore, / M is the filling 
fraction along the /i-direction, defined in Eq. |^. In this 
way, we can find the distribution of vortex segments n(r), 
by calculating the counterclockwise line integral of the 
currents around every plaquette in the system. Hence, 
the distribution of gauge-invariant current also gives in- 
formation, essentially by a duality transformation, about 
the FLL structure function. 

To perform a finite size scaling analysis, we employ 
the following procedure to find the current pattern in 
the ground state. Given a density of flux lines / = 1/32, 
we design an 8 x 8 vortex lattice unit cell, not to be con- 
fused with the unit cell of the numerical lattice. This 
vortex lattice unit cell has two vortices, Fig. [I]. The cur- 
rent pattern (Fig. [j]) is found by requiring the currents 
to obey Eq. ^ on every link and Eq. |] on every plaquette 
throughout the vortex lattice unit cell. It is possible to 
reduce the number of unknown currents by requiring the 
current pattern to have the same symmetry as the ground 
state vortex lattice. Periodic boundary conditions at the 
boundaries of the vortex lattice unit cell are used. 
By repeating the vortex lattice unit cell, we can design 
the current pattern of all systems with size m8xn8xl, 
(m,n,l) being positive integers. 

The flux-line density along the //-direction is defined 
as ffj,, and is given by 



B. The Ground State 



The current corresponding to the gauge invariance 
phase differences 

j M {f f ) = 0{i f +e lt )-0(f)-A pi {f f ), 

is defined on the directed link between site r and site 
7? + I -71- ! 71 ")- This current obeys two condi- 

tions: (1) There are no net current sinks or sources in 
the ground state 



V„j„(r)=0. 



(3) 



(2) The counterclockwise line integral of the currents 
around any plaquette of the numerical lattice with a di- 
rected surface normal in the /i-direction at site r must 
always be 



J2jAr) = 2n(n^r)-U). 



(4) 
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In the ground state, a uniform magnetic induction along 
the crystal c-axis B = Bz gives a periodic structure of 
straight flux lines aligned with B with hexagonal symme- 
try on a continuum substrate ab-plane, the well known 
hexagonal Abrikosov vortex lattice. In terms of the above 
densities / M this is expressed as 
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/, f* = f v = 0. 
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FIG. |l| a) The ground state flux-line configuration for f=l/32 
and system sizes L=8, 16,24. b) A cubic unit cell with two 
elementary vortex segments penetrating two plaquettes of the 
unit cell, c) The ground state current pattern in the ab-plane 
in units of 27r/384 for the 8x8 vortex lattice unit cell. The 
arrows indicate the direction of the currents on each link. 
The current pattern in this ground state is complicated, but 
nevertheless exhibits a high degree of symmetry. 

In our simulations, it is not possible to exactly load 
the hexagonal Abrikosov vortex-lattice onto our numer- 
ical mesh, which we have chosen to be square. This 
means that the underlying numerical mesh necessarily in- 
troduces a distortion of the hexagonal ground state. The 
numerical mesh represents a commensuration potential 
which acts as a perturbation on the ground state, and 
tends to 'freeze' the flux lines into a structure commen- 
surate with it. The flux lines will however tear them- 
selves off such a commensuration potential caused by the 
numerical lattice at a high enough temperature, which 
we denote a 'depinning' temperature Td p . Note that this 
depinning temperature has nothing to do with a real pin- 
ning potential, it is purely an artifact of the underlying 
numerical lattice. In the continuum limit, it would be 
zero. It is at present unclear to what extent the nu- 



merical mesh represents a singular perturbation on the 
continuum limit in a 3D system. The well known re- 
sults of Nelson and HalperincJ and YoungcJ concerning 
the effects of periodic commensuration potentials in 2D, 
indicates that if the filling fraction is small enough in 2D, 
Tdp will be smaller than any other relevant energy scale 
in the problem. This has also been nicely confirmed in a 
number of recent simulations on 2D systemso. 
Note that, by using a square numerical mesh, we coun- 
teract the disadvantage of the distortion of the hexagonal 
lattice by a reduction in the strong commensuration ef- 
fects we would have encountered if we had chosen a trian- 
gular numerical mesh, which admits an exact hexagonal 
lattice ground state. 

By using low enough filling fractions f z , it may be hoped 
even in 3D to achieve a satisfactory approximation to the 
continuum limit. That is, we hope that the depinning 
temperature Td p , which appears purely as an artifact of 
defining the model on a numerical lattice, drops below 
all other relevant temperatures in the problem, including 
the putative melting temperature T m of the FLL. 
That such a thermal depinning from the numerical lat- 
tice can actually be achieved in higher dimensions than 
2D is by no means clear, since commensuration effects 
are much more pronounced in 3D than in 2D. In fact, 
in the thermodynamic limit, thermal depinning from the 
numerical lattice strictly speaking cannot happen in 3D. 
Therefore, in order to mimic the statistical mechanics 
of the FLL defined on a continuum substrate in 3D, one 
carefully has to choose sample-geometry in simulations as 
follows: It is crucial to have a sample-geometry where the 
distance between flux lines is tailored to suit the thick- 
ness of the sample in such a way that the commensura- 
tion potential along the direction of the flux lines does 
not cause pinning to the numerical lattice at the tem- 
peratures of interest. This implies that simulations must 
be carried out on relatively flat slabs. Simulations on 
'tower'-shaped slaballl are probably not able to capture 
the continuum limit. In our simulations, we have checked 
that in fact such a depinning transition from the numer- 
ical mesh takes place at a temperature which is be- 
low the temperatures of primary interest. Hence, in our 
simulations the continuum limit ought to be adequately 
mimicked. Note that for filling fractions / = 1/8 and 
/ = 1/15, the continuum limit is not mimicked satis- 
factorily. In these cases, it is clear that pinning to the 
numerical mesh strongly influences the results. 
In this paper, we consider mainly the case of / = 
and / = 1/32, while some results are also obtained for 
/ = 1/72. 



C. The Monte Carlo Simulation 

The statistical mechanics of the Villain-model of 
an extreme type-II superconductor is investigated 
by employing the following Monte Carlo procedure 
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on numerical cubic lattices with linear dimensions 
L e [8,16,24,32,40,48,64,80,96]. 

Identical sets of current patterns are loaded onto each 
layer of the numerical lattice. For the filling fraction 
/ = 1/32, this current pattern is illustrated in Fig. [l]. We 
update the system, heating the system from the ground 
state consisting of fL 2 straight field-induced flux lines. 
A site of the numerical lattice is chosen randomly, and 
an attempt is made to change the phase on that site with 
a random amount A9 S [— ir, 7r). The phase change is ac- 
cepted or rejected according to the standard Metropolis 
algorithm. 

If the accepted phase change causes the current on a link 
jn(r) to fall outside the range j^{f) £ [— 7r, n) we add an 
amount ±2tt to the current, such that jfj,(f) is brought 
back into the primary interval j^(r) G [— tt, it). An im- 
portant point is that this operation can only generate a 
closed unit vortex loop around the link where the current 
is changed, thereby conserving the net induction of the 
system. No net vorticity is ever introduced by the pro- 
cedure, and the procedure also guarantees that no flux 
line can start or end within the sample. It is also impor- 
tant to note that the Monte Carlo procedure described 
above satisfies detailed balance. Hence, the entire phase 
space of the Villain-model is guaranteed to be exhausted 
provided the simulations are run for a long enough time. 
Another point is that the above procedure for limiting the 
currents to the primary interval also limits the number 
of vortex segments penetrating a plaquette to at most 
one per plaquette. In this way, the Villain-model dif- 
fers from the Lattice London Model, where the number 
of vortices penetrating each plaquette can be arbitrary. 
This difference however only becomes important in the 
high-temperature regime. 

The Monte-Carlo procedure really updates the gauge- 
invariant phase-differences, or currents j^(f). The simu- 
lations are therefore carried out in a manifestly gauge- 
invariant manner. One Monte Carlo sweep consists 
of L 3 attempts to change the phases 8(f) on L 3 ran- 
domly chosen sites throughout the lattice. Thus, by such 
a move we simultaneously change the gauge-invariant 
phase-differences on the six links associated with the 
relevant lattice point. Each data point for the quan- 
tities we consider is obtained after discarding the first 
100000 (30000) sweeps for equilibration. The subsequent 
400000 (70000) sweeps are used to obtain averages. The 
numbers in parentheses represent the number of sweeps 
we have for the case of T = 1. To ensure that measure- 
ments are independent of each other, we do one measure- 
ment per 100 sweeps. 



D. Anisotropy and finite size-effects 

For B — and isotropic couplings, the LSM, the XY 
Model, and the Villain Model all have one single phase 
transition at /cbT c /Jq ~ 3.0. The transition is charac- 



terized by stiffness in the phase of the superconducting 
order parameter being lost across the system in all direc- 
tions, due to a blowout of thermally excited closed vortex 
loops. 

For the anisotropic case, the bare coupling between 
planes J± in the Villain Model is smaller than the in- 
plane coupling J|| . Thus, in the very anisotropic case, 
the excitation energy of a unit vortex loop parallel to the 
ab-plane is much smaller than the excitation energy for 
a unit vortex-loop containing segments perpendicular to 
the ab-planeEl. One would naively then expect that ther- 
mal excitation of vortex-loops parallel to the ab-plancs 
would occur at correspondingly lower temperatures than 
those for which vertical loops appear. This is true for 
unit vortex loops, but such loops are unimportant for 
critical behavior. From the point of view of considering 
B = critical phenomena, the important issue is how 
the anisotropy affects large vortex loops, including vortex 
loops of order the system size. This is an issue to which 
we now briefly turn. It is convenient to carry out this 
discussion in terms of the phases of the superconducting 
order parameter, rather than in terms of vorticities. 
When the temperature approaches the Kosterlitz- 
Thouless (KT) temperature from above in a quasi-2D 
system, the phase-coherence length gradually grows. In 
a strictly 2D system, it would diverge precisely at the 
KT-transition. However, as long as a small coupling be- 
tween planes exists, no matter how small, then as the KT- 
transition is approached from above, increasingly larger 
domains of correlated phases are coupled together by the 
inter-plane coupling. This strongl^ptenormalizes the bare 
inter-plane coupling constant Jjj 15 H 22 l Hence, the system 
is isotropized close to, but above, the KT-transition, and 
the transition retains a 3D character. Thus, even an ex- 
tremely anisotropic system exhibits, precisely as in the 
isotropic case, one single 3D phase transition. No de- 
coupling transition as proposed in Ref.Eij exists in zero 
magnetic field. 

In simulations on finite systems, care must be taken to 
ensure that this physics is captured correctly. When 
the anisotropy of the system is increased such that the 
bare interplane coupling is reduced, one must make sure 
that the dimensions of the system in the directions par- 
allel to the ab-planes are large enough to allow the re- 
quired renormalization of the coupling constants to run 
its course without being cut off prematurely by the sys- 
tem size. Equivalently, one may rescale the size of the 
system in the z-direction, by the factor 1/r. Thus, the 
transverse system size must be tailored to the anisotropy 
of the system in such a way that critical behavior we 
study takes place at a lower temperature than the en- 
ergy scale ksT* set by 

Here, J± is the bare inter-plane coupling in the Villain- 
model, and £ is the coherence length of the phase of the 
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superconducting order parameter at the relevant temper- 
ature. If the system is too small in the transverse direc- 
tion, the renormalization of the bare coupling is cut off 
by the system size 

k B T* = L 2 J x . 

Hence, for a given L, we may choose an anisotropy such 
that J± is so small that T* becomes smaller than the 
actual temperature of the 3D critical phenomenon of in- 
terest, namely the vortex-loop blowout. We would then 
observe a decoupling of planes due to a proliferation of 
vortex-loops in the ab-plane which would be unphysical. 
This finite-size effect limits the anisotropies we can study 
consistently, at least in zero magnetic field. In zero field, 
we find ourselves limited to anisotropies of T < 4. 
A final technical point is that, although a finite magnetic 
field a priori allows larger anisotropies to be studied, the 
Villain-potential itself becomes virtually featureless as a 
function of its arguments Xi Eq. |2| for large anisotropies. 
Hence, simulations become impossible to perform mean- 
ingfully. 



in the continuum limit, T x and T y should vanish at all 
temperaturesE3, since a current in the ab-plane would 
exert a Lorentz-force on the unpinned flux lines, mov- 
ing them and thus dissipate energy. When the model is 
discretized by introducing the numerical lattice, a finite 
energy for moving them in a direction perpendicular to 
the z-dircction is also introduced. The existence of a 
smallest energy to required move flux-lines acts as an ar- 
tificial pinning potential on the FLL, causing T x (T y = T x ) 
to be finite. Thus, the FLL 'depins' from the underlying 
discrete lattice at a finite temperature T x = T y > 0. In 
the continuum limit, as long as no physical pinning of 
the flux lines is present, we would have T x = T v = 0, 
and the flux lines are unpinned at all temperatures. To 
ensure that the above artificially introduced pinning po- 
tential caused by the numerical lattice does not affect 
the FLL melting transition and any other genuine phase- 
transition we might want to consider, we should consider 
systems with T m significantly higher than T x . The way 
to achieve this is to consider low enougi .filling fractions 
/ = f z of flux lines. Several authorsEJo have in fact 
found that T x decreases for decreasing flux lines density 
f, and falls below T m for / < f c ~ 1 /32. 



E. The helicity modulus 



To probe the global phase coherence in the BCS su- 
perconducting order parameter across the entire system, 
we consider the helicity modulus T M , defined as the sec- 
ond derivative of the free energy with respect to a phase 
twist in the /i-dircctionE£l. It basically measures the stiff- 
ness of the system to a twist in the phase of the order 
parameter. In the anisotropic case we have the following 
generalization of previously obtained expressions for the 
helicity modulus for isotropic superconductorstll 




For the details on the derivation of this expression, and 
a corresponding one in terms of vorticities, see Appen- 
dices A and B. Here, V,', and V" are the first and second 
derivatives, respectively, of the anisotropic Villain poten- 
tial V/j, defined in Eq. |2|. For temperatures T < T M such 
that T M > 0, there is phase coherence across the entire 
system in the /i-direction. Hence, the system can sustain 
a supercurrent in the ^-direction. At the temperature 
T = T M and above, such that T M = 0, phase-coherence 
along the ^-direction is lost. Hence, the vanishing of 
T„ signals the superconducting normal metal transition 
in the //-direction with the transition temperature T^. 
For the case of finite magnetic induction along z-axis 



F. Vortex-Loop Distribution 

As mentioned in the Introduction, in the LSM, Villain- 
model, and Lattice London Model, the zero-field normal 
metal-superconductor transition corresponds to a vortex- 
loop blowout analogous to what has long been suggested 
to occur in the neutral superfluid He*. To study the 
blowout of closed vortex loop in extreme type-II super- 
conductors, we consider the quantity D(p), which we de- 
note the vortex-loop distribution function, and which is 
given as statistical average of the total number of closed 
non-field induced vortex loops with a given perimeter p, 
in our case normalized by the volume of the systems we 
consider. 

The following procedure is employed to compute D(p). 
We start from an arbitrarily chosen unit cell containing 
at least one outgoing vortex segment penetrating a pla- 
quette of that unit cell. We then follow the direction 
of this vortex segment into the neighboring unit cell. If 
there is more than one vortex segment leaving the unit 
cell, one of them is chosen randomly. We continue trac- 
ing the path of vortex segments until the path closes up 
on it self. When the path has closed upon itself, we mea- 
sure the length I of the path, as well as the net vorticity 
v z along z-axis of the path. Also, we remove the vor- 
tex segments along the path to prevent double counting 
of paths. Because of the periodic boundary conditions, 
such closed paths of vortex segments can either belong 
to a field-induced flux line or a non field induced closed 
vortex loop. 

In the quantity D(p), we do not include the closed paths 
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associated with field-induced flux lines that close on 
themselves merely due to periodic boundary conditions 
in the z-direction. Field-induced flux-lines are character- 
ized by a net vorticity in the z-direction, v z ^ 0. For the 
purposes of studying the loop-transition, we are exclu- 
sively interested in closed paths associated with non-field 
induced vortex loops that physically close on themselves 
regardless of boundary conditions. Hence, the relevant 
closed paths of vortex segments are closed vortex loops 
with perimeter p = I which have net vorticity = in 
all directions. 

The procedure for tracing out the relevant closed loops is 
repeated until all non-field vortex segments in the system 
have been counted. For B ^ 0, this procedure uniquely 
separates thermally excited closed vortex loops from the 
field induced flux lines. 

In the low temperature regime D(p) depends on the exci- 
tation energy E(p) ~ ep of the -saartex loops with perime- 
ter p, with a Boltzmann-factortZI, 

D(p)~exp(— )~exp(— p-), 
k B l L Q 

where £ is a constant representing a line tension, and 
Lq ~ kgT/e is a typical perimeter of closed vortex- loops 
present at a given temperature in the low-temperature 
regime. As we will see below, such low-temperature "con- 
fined" vortex-loops may be coarse-grained away and are 
unimportant for the statistical mechanics of the mixed 
state of an extreme type-II superconductor. In this ther- 
mally activated regime, large vortex loops are exponen- 
tially suppressed. On the other hand, at the critical 
point, vortex loops with all perimeters are present. This 
leads to an algebraic decay of the loop-distribution func- 
tion versus loop-perimeter at the critical point, 
D(p) ~ p~ a , where a is an exponent not to be confused 
with the critical exponent of the specific heat. Hence, 
monitoring the temperature where D(p) changes its char- 
acteristic behavior from exponential decay to algebraic 
decay, is a way of determining the vortex-loop unbinding 
temperature. 

If we assume that the vortex-loop-distribution function 
scales with the vortex-loop perimeter as some power-law, 
and furthermore assume that the perimeter scales with 
the vortex-loop radius r, then we have 

D(r) ~ r- a . 

We now use a critical scaling analysis to determine a 
in our case. The assumption is that the loop-transition 
in zero field represents a critical point. If we can then 
fit numerically obtained exponents at the putative criti- 
cal point to the scaling results, this would provide further 
support for the assertion that the loop-blowout is respon- 
sible for destroying superconductivity in extreme type-II 
superconductors. 

The Villain-model is dual to a 3D Coulomb-gas, and the 
vortex-loops are analogous to the vortex-antivortex pairs 
of the 2D Coulomb-gas. We may determine the contri- 
bution to the dielectric constant of the 3D Coulomb-gas 



that such loops give. Since the dielectric constant may 
be related to the inverse superfluid stiffness, whose scal- 
ing dimension is well known on general grounds, we may 
determine a at the critical point. 

The "dipole-moment" P(r) of a vortex-loop scales with 
r d-2 , r ^ w h ere the charge ~ r d ~ 2 , and the dipole vec- 
tor ~ r. The contribution to the dielectric constant, or 
the electric susceptibility, coming fpp thermally induced 
loops of size between r is given byE3 

d 

Xe(r) ~ e(r) ~ < r d ~ x cos(</>) >\ E =o, 

where the average should be a thermal average with the 
Boltzmann- factor 

exp[-U(r)/k B T] ~ D(r) exp [r d_1 E cos(»/fc B T]. 

Here E is an electric field polarizing the medium via the 
"charge-loops" of the 3D Coulomb-gas, and (j> is the angle 
between the orientation of loops and the applied electric 
field polarizing the medium. We find 

On the other hand, in the superconductor, the super- 
fluid stiffness p s is given by the transverse susceptibility 
X ~ G ~ p s , where G is given by the order-parameter 
Green's function, G ~ r 2 ^ d ^ r) , and where r\ is the anoma- 
lous scaling dimension of the Green's function appearing 
due to critical fluctuations, for the 3-DXY-model we have 
r\ = 0.033(4). Now we use the fact that e ~ pj 1 , to find 

2{d- 1) - a = d + ri -2, 
a = d — rj. 

In our case, we may evaluate the loop-distribution func- 
tion at the anomalous peak in the specific heat, and 
fit the result to a power law with the exponent 3, ob- 
tained by assuming critical scaling. We will see below 
that the fit is excellent, lending further support to the 
assertion that in zero field, the superconductor-normal 
metal phase-transition in an extreme type-II supercon- 
ductor is due to a vortex-loop blowout transition. 
It is also interesting to note that another way of estimat- 
ing the relevance of closed vortex-loops, is to see whether 
they can be coarse-grained away or not. A rough cri- 
terion for coarse-graining them away, would be that the 
loop-distribution integrated over the volume of the sys- 
tem should be finite, i.e. 

J d d r D(r) ~ J dr r^ 1 -" ~ A d ~ a . 

Thus, loops may be coarse-grained away provided that 
a > d. In the low-temperature regime, we have seen 
that the loops are even exponentially suppressed, and 
certainly satisfy this criterion. The marginal case a — d 
gives an integrated distribution ~ In A. However, at the 
critical point, a = d — rj is less than the required value 
for coarse-graining, furthermore critical fluctuations as 
manifested by a non-zero positive value of the exponent 
rj, will increase the relevance of vortex-loops, as expected. 
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G. Specific heat 

In addition to indirect measurements of the latent 
heat of melting, such as measured by local magnetization 
measurements on BSCCOca, direct calorimetric measure- 
ments of the specific heat are also useful for estimating 
the latent heat of the melting transition of the FLL, or 
any other phase-transition the vortex system might suf- 
fer. Such measurements ar.©-|now availabled both in zero 
field and in finite field. InO, the specific heat anomaly 
in a twinned YBCO sample is measured systematically 
with varying magnetic field. It evolves smoothly from 
the zero field result as the magnetic field is increased. 
Moreover, the integrated anomaly appears to be approx- 
imately constant as the magnetic field increases. 
This raises the question of what sort of phase-transition, 
if any, the specific heat anomaly in finite magnetic fields 
should be associated with. Due to its smooth evolution 
from the zero-field case, it appears rather unlikely that 
this sizeable anomaly has anything whatsoever to do with 
FLL melting. Rather, it seems to suggest that there are 
remnants of the zero-field transition, which we will de- 
scribe in detail below, at finite magnetic fields. 
We may investigate this issue, by calculating the specific 
heat of the Villain model, and correlate the specific heat 
anomalies with the temperature dependence of the struc- 
ture function S(K) as well as with the phase-stiffness T z . 
Thus we should in principle be able to decide whether or 
not the major features in the specific heat have anything 
to do with FLL melting or vortex- loop blowout, also at 
finite fields. In zero field, we will be able to precisely 
correlate, for all anisotropies, the anomaly in the specific 
heat with the vortex-loop blowout transition. In a finite 
field, the situation is considerably more complicated. We 
find three anomalous features in the specific heat. The 
major feature in the specific heat, the remnant of the 
zero-field transition, occurs at temperatures well above 
those where the structure function of the FLL, and the 
phase-stiffness across the sample along the field direction, 
vanishes. 

To calculate the specific heat per site C, we use the usual 
fluctuation formula 

C_ _ I <H 2 V > - < H v > 2 

~ 1? {k B Tf " 

The Villain model has a rather unusual property in that 
the Boltzmann factor appearing in the partition function, 
exp(— H/ksT) involves an explicitly temperature depen- 
dent Hamiltonian, Eq. (En. The usual fluctuation for- 
mula for the specific heat is valid strictly speaking only if 
H is temperature independent. Nevertheless, we will use 
the above standard expression for calculating the specific 
heat for convenience, and neglect the extra terms that 
should be included from the explicit temperature depen- 
dence of the Villain potential, Eq. g. We have checked the 
validity of this approximation by comparing the thus ob- 
tained specific heat per site C with the alternative stan- 
dard method of extracting the specific heat per site C 



from its basic definition in terms of the internal energy 
of the system 

C I d <H V > 
k^~U d(k B T) ■ 

We find that these two ways of calculating the specific 
heat differs only in the very high temperature regime, 
outside the temperature range of interest in this paper. 



H. The Structure Function 

To probe the FLL melting, we consider the structure 
function for n z vortex segments, i.e. vortex segments di- 
rected along the averaeeinduction. The structure func- 
tion S(k) is defined by00 

ci/jA _ <l E^(f)exp [ik-r\ | 2 > 

For our ground state with the flux lines density / = 1/32 
(Fig. |f|), the unit reciprocal lattice vectors for the FLL 
are 

j? 1 = 27r[i,~], K 2 = 2n[0,±\. 

In the FLL phase, S(k±,k z = 0) has (5-function 
Bragg peaks at k± = K(m,n) = mK\ + nK 2 
(m,n=0,±l,±2,±3,..). The vectors K are located within 
the first Brillouin-zone. When the FLL melts, the Bragg 
peaks are smeared out. The lowest temperature T where 
S(K, k z = 0) vanishes, thus defines the FLL melting tem- 
perature T m . For simplicity, we consider only the struc- 
ture function S(Q) = S{K 2l k z =0). 



I. Flux-line cutting and intersection 

The issue of flux-line cutting and flux-line and entan- 
glement, as well as the suggested possible resulting vortex 
states originating from the latter, such as the analogs of 
2D Bose superfluids and supersolids, and even topolog- 
ical vortex glass states, have been issues of considerable 
controversy over the last years. In particular, the effect 
of entanglement on the FLL melting transition and the 
statistical mechanics of the FLL, has received consider- 
able attention. 

In principle, flux-line entanglement could be responsible 
for the drop in T 2 we observe at the temperature T z . It 
is of interest to correlate the amount of "close vortex-line 
encounters" with the anomalies we obtain in the specific 
heat. This will allow us to at least tentatively decide 
whether or not flux lines start to intersect or cut at any 
of the temperatures T m , T z , or Tb c i- 
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Intersection and cutting of flux lines tend to act as effi- 
cient modes of disentanglement of flux lines. It is unlikely 
that a flux-line liquid phase that suffers large amounts of 
thermally induced collisions between flux-line segments, 
which in turn strongly indicates that flux-line cutting 
takes place, can also sustain heavily entangled vortex 
configurations. Hence, if we can show that the amount 
of close vortex-segment encounters is substantial at the 
temperature where T z goes to zero, we may at least ten- 
tatively conclude that the drop in Y z is not due to en- 
tanglement of flux-lines. 

To study the intersection and or cutting of flux lines, we 
define a parameter p, which is a measure of the amount 
of flux-line cutting and intersection that takes place in 
the flux-line liquid. 



where Ni nt is the total number of unit cells having 4 or 
6 vortex segments penetrating their plaquette, and is 
the number of field induced flux-lines in the system. Note 
that with this definition, we only consider p for the case 
of finite fields. This quantity gives us some intuition on 
how well-defined we may expect individual field-induced 
flux lines to be in the liquid phase. 

In the Villain Model, each plaquette can carry at most 
one vortex segment. This is because the vorticities 
are defined in terms of gauge-invariant phase-differences, 
defined on the interval jn(r) G [— it, tt). The dis- 
tribution of vortex segments n(r) has no divergence 
(E^x.^z^^mW = °)- Thus, a single unit cell can 
only carry 0, 2, 4, or 6 vortex elements. If a unit has 
more than 2 plaquettes that are pierced by vortex seg- 
ments, then that corresponds to one of the three following 
situations: 

(1) an intersection between two different flux lines, or 

(2) an intersection between a flux line and a closed vortex 
loop, or 

(3) an intersection between two different closed vortex 
loops. 

Inside a unit cell, it is impossible to decide which possi- 
bility is realized. In principle, p as defined above counts 
all of these possibilities, while flux-line intersection and 
flux-line cutting correspond to case 1). 

III. RESULTS, ZERO MAGNETIC FIELD 

We now present the results of our Monte-Carlo simu- 
lations, and consider first the case of zero magnetic in- 
duction B = 0. From now on we measure the specific 
heat per site C in units of fcs, the helicity moduli T x in 
units of Jo, T z in units of Jq/T 2 , and the temperature T 
in units of Jo/fc_B. 

Two values of T = A C /A Q , the anisotropy parameter, are 
considered: T = 1 and T = 3. We will present results for 
the quantities relevant to the zero-field case, namely the 



helicity modulus, the vortex-loop distribution function, 
and the specific heat. An important point is that we need 
to, in the zero-field case, to be able to correlate the tem- 
perature at which the helicity modulus T z vanishes, with 
the temperature at which the loop-distribution function 
D(p) qualitatively changes behavior from an exponen- 
tial dependence on the loop-perimeters to algebraic de- 
pendence on the loop-perimeters. Moreover, both these 
features must be correlated with the temperature where 
an anomaly in the specific heat is found, as discussed 
in Section II. D. For arbitrary anisotropy T, the system 
should only have one single phase-transition, the normal 
metal-superconductor transition. 



A. Helicity modulus 

The results for the helicity modulus T M , Eq. |^, in zero 
magnetic field and in the isotropic case T = 1, are shown 
in Fig. |^. We have confirmed that all moduli T X1 T V ,T z , 
are equal in this case, and therefore only exhibit T z . The 
results are shown for the system sizes L — 8, 32, 64, 96. 
Note that the drop in T z becomes sharper as the system 
size increases, and the value of T where T z appears to 
vanish, becomes smaller. When A — * oo, T z , which is the 
stiffness of the phase of the superconducting order param- 
eter to a twist, is proportional to the superfluid density 
p s . (For a finite A this identification no longer holds, as 
emphasized inE3) . We have found from the numerics that 
T z ~ \T — T c | 2 / 3 , consistent with the Josephson scaling 
relation, p s ~ \T — T c \ u , where v is the correlation length 
exponent v = 2/d = 2/3. 




2 2.3 2.6 2.9 T 3.2 



FIG. [2| Helicity modulus Y z and the specific heat per site 
C versus temperature, for B = 0, F=l (isotropic) and system 
sizes L=8,32,64,96. 
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Due to the identification Y z ~ p Sy we conclude that 
the vanishing of the helicity modulus corresponds to the 
normal metal-superconductor transition. The transition 
occurs at the temperature T = 3.0 in units of Jo/fee- 
The results for the helicity modulus in the anisotropic 
case r = 3 are shown in FigJ^. The situation at first 
glance appears considerably more complicated than in 
the isotropic case, despite our expectations that the 
physics basically should be the same as in the isotropic 
case, cf. our discussion in Section II. D. A striking fea- 
ture is that for r — 3, the helicity modulus Y z appears 
to vanish distinctly below the temperatures at which the 
helicity moduli , T y vanish. (T^ and T y turn out iden- 
tical in all our simulations, and we therefore only exhibit 
T x ). Note however that there is an important finite- 
size effect in the results: As L increases, vanishes 
at progressively lower temperatures while T z vanishes at 
progressively higher temperatures. As L increases, Y z 
and T x appear to approach zero at the same tempera- 
ture. Due to the limitations in available system sizes, we 
have not been able to perform the simulations at higher 
anisotropies than r = 3 for the zero-field case. At lower 
anisotropies 1 < T < 3, the same finite-size effect as de- 
scribed above is seen. For lower values of the anisotropy it 
is also more obvious that the two temperatures at which 
T z and T x vanish, approach each other with increasing 
system size L. For the anisotropic case T = 3, the tran- 
sition occurs at T = 1.57 in units of Jq/^s- 
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FIG. |3|. Helicity moduli T x , T z and the specific heat per 
site C versus temperature, for B = 0, T=3 and system sizes 
L=8,32,64,96. For increasing system size L, T x approaches 
T c — 1.57 from above, and T z increases approaching T c from 
below. The layered system has only one phase transition. 



B. Loop distribution 

In order investigate the excitations responsible for de- 
stroying the superconducting phase-coherence and the 
superfluid stiffness as evidenced by our results for T z , 
we probe the amount of closed vortex-loops that are ther- 
mally excited in the superconductor model at the temper- 
ature where Y z vanishes. We first discuss the isotropic 
case r = 1. 

The results for D(p) are shown in Fig. || for the largest 
system we have considered, L = 96. The figure shows 
D(p) as a function of the loop perimeters p, for various 
temperatures in the range T £ [2.3,3.3]. Recall from 
above that in the isotropic case, the helicity modulus 
vanished at T k 3.0. 
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FIG. i. The distribution function D(p)/L 3 of thermally ex- 
cited, non-field induced, closed vortex loops per site as a 
function of perimeter p, for B = 0, T=l, L=96 and various 
temperatures. It has been normalized this way to facilitate 
comparison between different system sizes. For T < T c ~ 3.0, 
D(p) is best fit to an exponential decay. For T — T c , D(p) 
is fits an algebraic decay D(p) ~ p~ a with exponent a ~ 3 
excellently, indicating an Onsager loop-transition at T = T c . 
The inset shows D(p) on a log-log plot. The slope of the 

straight line, obtained at T = T c , is ' 3. At T < T c 

the curves show a marked downward curvature, indicating a 
faster-than-algebraic decay of D(p). This point is discussed 
in section II. F. 

From Fig. ^, particularly from the inset of this figure, 
we observe a qualitative change in D(p) precisely at the 
temperature T = 3.0. The inset shows the distribution 
function on a log-log plot, and it is seen that the decay 
is faster-than-algebraic for T < T c while it is a precise 
power law with exponent a in good agreement with a 
scaling analusis assuming that the vortex-loop blowout 
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is a critical point. We attribute the slight deviation in 
the exponent a between the simulations and the theory 
as due to the presence of vortex loops of more compli- 
cated shapes, such that the circumference of the loops 
do not all scale with their diameter. We have tentatively 
suggested an exponential decay in the low-temperature 
phase, but this is not unambiguous. However, our main 
point is that for T < T c the decay cannot be a power law, 
while at T = T c , the power law we find is precisely the 
same as the one we predict analytically assuming that 
the vortex- loop blowout is in fact a critical point. 




FIG. ||. The distribution function D(p)/L 3 of thermally 
excited closed vortex loops per site as a function of perimeter 
p, for B = 0, r=3, L — 96 and various temperatures. For 
T < T c ~ 1.6, D(p) is best fit to an exponential decay. For 
T = T c , D(p) is best fit to an algebraic decay, D(p) ~ p~ a , 
with a ~ 3 consistent with critical scaling of section II. F. The 
inset shows the same figure for L=16. Note that for L=16, 
D(p) incorrectly shows algebraic decay for T = 1.3 < T c . 
This finite-size effect is discussed in Section II. D. 

Hence, based on the above, we conclude that at 
T = 3.0, a sharp phase-transition occurs from a low- 
temperature phase where closed vortex-loops are con- 
fined to some typical size Lq(T), to a high-temperature 
phase where closed vortex loops of all sizes up to and 
including the system size, exists. Thus, in the isotropic 
case we have been able to precisely correlate the drop in 
T z with a vortex-loop blowout, and from the previous 
paragraph we must identify this as the superconductor- 
normal metal transition. 

The loop-distribution function for the anisotropic case 
r = 3 is shown in Fig. ||. Due to the drop in the crit- 
ical temperature of the system, we now show D(p) as a 
function of p for temperatures in the range T G [0.7, 2.1]. 



Again, we observe a qualitative change in the behavior of 
D(p) from exponential decay to algebraic decay, at a tem- 
perature T w 1.6, which correlates almost perfectly with 
the temperature T — 1.57 at which the helicity moduli 
T z and T x vanish for T = 3. If we fit D(p) ~ p~ a at this 
temperature, we again find the exponent a = 3, as in the 
isotropic case. 

So far, our expectations based on the discussion in Sec- 
tion II. D are borne out. To illustrate the point further, 
in the inset of Fig. ^, we show the distribution function 
D(p) for r = 3 and the same range of temperatures, for 
the smaller system L = 16. The important difference 
between these two cases is that for L = 16, algebraic 
decay of D(p) appears to persist down to lower tempera- 
tures than for L — 96. For larger systems the vortex-loop 
blowout is suppressed due to the fact that the interplane 
coupling is allowed to renormalize further without being 
cut off by a small system size. Hence, what appeared to 
be a separate vortex-loop blowout at a low temperature 
T = 1.1 for L = 16, has been pushed up to the correct 
temperature T « 1.6 in the larger system L = 96, as 
discussed in Section II. D. 



C. Specific heat 



We now present our results for the specific heat in zero 
magnetic field for the system sizes L = 8, 32, 64, 96, and 
consider first the isotropic case T = 1. 
The specific heat in the isotropic case is shown in Fig. 
||. The anomaly in the specific heat clearly correlates 
with the temperature where the supcrfluid stiffness van- 
ishes, which in turn correlates precisely with the temper- 
ature where the vortex-loop blowout is observed. As the 
system size increases the anomaly clearly also becomes 
sharper. We have shown that the peak in the specific 
varies in very good agreement with ln(i) as the system 
size is increased, further indicating a genuine thermody- 
namic phase transition. The shape of the specific heat 
curve has a typical Ay-behavior for this extreme type- 
II superconductor (A — > oo). This agrees with previous 
results found by Dasgupta and Halperin for the lattice su- 
perconductor modelE3. (Note that this contrasts sharply 
with the specific heat results indicating an inverted XY- 
transition found by the same authors for finite, small A in 
the isotropic case). Qualitatively, our results also agree 
well with the [Specific heat measurements on YBCO of 
Schilling et alEB. 

The specific heat for the anisotropic case T = 3 is shown 
in Fig. U for the system sizes L = 8, 32, 64, 96. The situa- 
tion again at first glance appears to be more complicated 
than the isotropic case. However, as in the isotropic case, 
the peak in the specific heat anomaly increases with sys- 
tem size. Moreover, the temperature of the peak in the 
specific heat is reduced as the system size is increased, 
approaching the temperature T = 1.57 at which both T 2 
and vanish, and where the vortex-loop blowout ap- 
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pears to take place. 

Hence, our simulations of the above three quantities T z , 
D(p), and specific heat correctly capture the physics that 
even in the very anisotropic case, only one single phase 
transition occurs in the lattice superconductor model in 
zero magnetic field. The fact that our simulations do not 
reveal an artificial zero-field decoupling transition due to 
finite-size effects, makes us confident that our simulations 
are able to capture the subtle zero-field physics correctly 
in the anisotropic case. This is a necessary prerequisite 
for being able to extract meaningful results from our fi- 
nite field-simulations, to which we now turn. 



IV. RESULTS, FINITE MAGNETIC FIELD 



Next, we present results for finite magnetic induction. 
We will consider the filling fraction / = 1 /32 correspond- 
ing to 32 ab-plane plaquettes per field-induced flux line 
in the ground state, depicted in Fig. |l|. Again, we are pri- 
marily interested in correlating temperatures where the 
helicity moduli vanish with temperatures where anoma- 
lies in specific heat occur, as well as with temperatures 
where we see a qualitative change in the distribution of 
non-field induced closed vortex-loops. It is also of inter- 
est to correlate these phenomena with the melting of the 
FLL as evidenced by a drop in the structure function at 
low-order Bragg-peaks. Furthermore, due to the presence 
of field- induced flux lines, it is also of interest to monitor 
the amount of flux-line cutting occuring in the system. 
This has bearing on the amount of entanglement of flux 
lines the molten vortex phase can sustain. For finite-fields 
we present results for helicity moduli, vortex-loop dis- 
tribution, specific heat, structure function and flux-line 
cutting for the two values of the anisotropy parameter, 
r = 1 and T = 3. 



A. Structure function 



The results for the structure function S(Q) are shown 
in the top panel of Fig. || for the isotropic case T = 1, 
for the reciprocal lattice vector Q = (K2, k z = 0), K% = 
2tt[0, 1/4], and for the system sizes L = 48, 64, 80, 96. In 
the top panel of Fig. fj], results for V = 3 and for various 
system sizes L = 32, 48 are also shown. The structure 
function exhibits a sharp drop at T = T m , well below the 
temperature where the phase-coherence in the direction 
of the magnetic field vanishes. (We discuss the helicity 
modulus T z in the next section.) Note that the structure 
function vanishes essentially at the same temperature as 
the temperature where the helicity modulus T x vanishes. 



L = 16 □ 




FIG. |6j. Top panel: Helicity moduli Y x and Y z and struc- 
ture factor S(Q) versus temperature, for / = 1/32, T=l and 
L — 16, 32, 48, 64, 96. Note how for increasing system size, T m 
increases slightly, and T z decreases markedly. Lower panel: 
The specific heat per site C versus temperature for / = 1/32, 
r=l and L = 16, 32, 48, 64, 96. 

What this indicates is that the filling fraction we have 
used in our simulations, / = 1/32, is not sufficiently small 
to study the melting of the FLL in the continuum limit. 
The FLL has not thermally 'depinned' from the numer- 
ical mesh at the melting transition. Therefore, our es- 
timate for the FLL melting temperature as obtained in 
these simulations is too high. However, from our ear- 
lier wpdc on the moderately anisotropic Lattice London 
modefcy, we know that / = 1/48 suffices to produce a 
thermal 'depinning' temperature of the FLL off of the 
numerical mesh at a temperature distinctly below the 
observed melting temperature of the FLL. Therefore, we 
expect that the present estimate for T m should be quite 
good; only a minor reduction of the filling fraction below 
/ = 1/32 is expected to suppress the 'depinning' temper- 
ature T x below the melting temperature T rn . It is also 
possible that the commensuration effect might tend to 
overestimate a first order character of the FLL melting 
transition, should such a result be found. 
The reduction in the structure function below the melt- 
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ing temperature is due to the Debye- Waller factor 

S(Q, T) = S{Q, 0) exp(-G 2 < u 2 >) (8) 

For a triangular lattice, we find right below the melting 
transition 

S(Q,T) = S(Q,0) exp(-^ci) (9) 

where cl is the Lindemann-ratio. In our simulation, we 
find that the DW-factor is 0.6 right below the melting 
transition, and hence cl = 0.24. Essentially the same 
result is found for the isotropic case, and is in reasonable 
agreement with the value cl = 0.4 used in the best calcu- 
lation so far to estimate the position of the FLL-melting 
line in BSC CO and YBCO by employing the simple 
Lindemann-criterion in conjunction with the highly non- 
trivial fluctuation propagator found from anisotropic and 
non-local elastic theory of the flux-line latticeEil. 



B. Helicity modulus 

The results for the helicity moduli T z and T x are 
shown for T — 1 and system sizes L — 16, 32, 48, 64, 80, 96 
in the top panel of Fig. |^. Likewise, similar results for 
r = 3 are shown in the top panel of Fig. for system 
sizes L = 16,32,48. Note how the temperature T z ap- 
pears to decrease monotonically with system size. 
An important issue is how T z will continue to vary when 
the system size is increased indefinitely. Since there is no 
obvious sign of saturation in T z as L increases, it could 
conceivably continue to decrease in the liquid phase, un- 
til it reaches T m . Does this in fact happen, or is T z > Tm 
in the thermodynamic and continuum limit? 
To answer this question, we have performed simulations 
on systems with filling fractions / = 1/72, and computed 
T m and T z as a function of L. The point about going to 
lower filling fractions is that we are approximating the 
continuum limit better. It is becoming increasingly clear 
from numerical simulations that when the xy-plane is 
discrctized in order to do the simulations one is intro- 
ducing a long time scale into the problem: There is a 
gap for moving vorticities from one unit cell to another. 
When the filling fraction is lowered, the continuum ap- 
proximation is better approximated, and the relaxation 
time introduced by discretization is lowered. Hence, for 
/ = 1/72 we are better able to equilibrate the system. 
(Incidentally, we believe the peason that no finite-size 
effect was seen in T z in RefO was precisely that the 
simulations were not run for a long enough time). The 
results of our simulations for this cases is shown in Figure 
|, where we show T Z (L) - T m (L) for V = 1 and T = 3, 
for the case / = 1/72. It is seen thatv T Z (L) - T m (L) 
decreases monotonically as a function of L. Is it possible 
that T z could drop arbitrarily far below T m when the 
anisotropy T is increased indefinitely? We believe that 



the answer for A = oo is no, for the following reason. 
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FIG. Top panel: Helicity moduli Y x and Y z and struc- 
ture factor S(Q) versus temperature, for / = 1/32, T=3 
and L = 16, 32, 48. Note how for increasing system size, T m 
increases slightly and T z decreases markedly. Lower panel: 
The specific heat per site C versus temperature for / = 1/32, 
r=3 and L — 16, 32, 48. Note the two separate specific heat 
per site anomalies at T x = T m and T z . The specific heat 
peak at T m becomes more prominent with increasing system 
size, while the peak at T z actually decreases with increasing 
system size. 

Recall our discussion in Section II. D, where it was 
shown that in zero field, no vortex-loop blowout can take 
place below the zero-field transition temperature, due 
to the renormalization of the interplane coupling. In fi- 
nite fields, this is quite different, since the superconduct- 
ing coherence length is limited by the magnetic length 
as soon as the flux-lines start to fluctuate appreciably, 
cutting off the renormalization of the interplane phase- 
coupling. Hence, we see that T z drops well below Tb C 2, 
in finite fields, provided the system is in the liquid-phase. 
In the flux-line lattice phase, the coherence length is no 
longer limited by the magnetic length, there is now phase- 
coherence throughout the sample (but of a more compli- 
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catcd form then in the zero-field case) . Hence, the renor- 
malization of the interplane coupling described in Section 
II. D again becomes active, suppressing the vortex-loop 
blowout. The flux- line lattice phase is therefore in some 
sense equivalent to the zero-field case with regards to a 
vortex-loop blowout, and a loss of phase-coherence along 
the field direction cannot take place within the lattice 
phase. Hence, T z cannot drop far below T m , at least not 
for the case A = oo. 



0.8 

0.7 

0.6 

2 0.5 
_E 

7 0.4 
0.3 



0.2 



f = 1/72 



r=i 
r=3 



lines is present. On the other hand, T z has little or no 
commensuration effects in it. 

The above result indicates that phase-coherence across 
the sample along the direction of the flux lines is lost 
at a finite temperature T z . Contrary to the zero- field 
case, there are several possible explanations for this loss 
of phase coherence when a finite field is present. One 
possibility is that a vortex-loop blowout causes the loss 
of phase-coherence at T = T z . Another explanation 
could be that since the melting temperature T m is smaller 
than T z , in FLL liquid phase flux lines become entan- 
gled thereby destroying the phase coherence in the su- 
perconducting order-parameter along the direction of the 
flux lines. A third explanation could be that transverse 
flux-line meanderings and flux-line cutting causes loss 
of phase-coherence without a resulting entanglement of 
flux-lines. We discuss these possibilities in turn. 



C. Loop distribution 
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FIG. |. The difference T Z (L) - T m (L) as a function of L 
for / = 1/72, T = 1 and Y = 3. Note that this quantity 
shows a monotonically (sublinear) decrease with L. Note also 
that the relative decrease is larger for F = 1 than for Y = 3. 
We attribute this to a slower relaxation in the anistropic case 
than in the isotropic, since the both T m and T z are smaller 
for T = 3 than for V = 1. 

In our previous work0 we emphasized the importance 
of closed loops, but suggested that the drop in T z within 
the liquid phase was a genuine phase-transition from a 
coherent vortex liquid characterized by T z 0, to an 
incoherent vortex-liquid characterized by T z = 0. We 
strongly believe our present simulations on much larger 
systems show that this in fact is not the case: T z — > T m 
as the system size increases. However, the final word on 
the issue, particularly for the isotropic case Y = 1, re- 
mains to be said. 

Note also that, compared to the zero-field case / = 0, 
T z vanishes at a considerably lower temperature when 
/ = 1/32. The helicity modulus vanishes at a finite 
temperature below T z . 

The fact that T x vanishes at a finite temperature is an 
artifact of our discretizing the ab-plane. In the contin- 
uum limit, this helicity modulus would be zero for any 
finite temperature when no physical pinning of the flux 



The results for the distribution of closed vortex-loops 
are shown in Fig. |9[ in the temperature range T 6 
[0.9,2.5], for the anisotropy Y = 3 and the system size 
L = 48. This temperature range encompasses the melt- 
ing of the FLL and the destruction of phase-coherence 
along the direction of the magnetic field. A feature 
which distinguishes the finite-field results for D(p) from 
the zero-field case, is that throughout the temperature 
range where the FLL melts, T m « 0.38 and where phase- 
coherence is lost, T z 0.65, the distribution function 
D(p) decays more rapidly than at the critical point in 
zero field. In the range T < 0.7 there is no obvious sign 
in D(p) of a cross-over to algebraic decay as a function of 
loop-perimeters. We conclude that the vanishing of Y z is 
not, in this case, associated with a vortex- loop blowout 
on all length scales of the system. In other word, the 
vanishing of T z is not due to a finite-field counterpart of 
the Onsager-loop blowout we found in zero field. We find 
a change in the behavior in D(p) from exponential decay 
to algebraic decay at a much larger temperature T w 2.8 
for r = 1 and T w 1.9 for Y = 3. The temperature 
range over which D(p) changes behavior is also consider- 
ably broader than in the zero-field case, indicating that 
the vortex-loop blowout transition which was found to be 
sharp at zero field, is replaced by a crossover. 
However, the interaction between closed vortex-loops and 
the flux-line lattice may be studied by considering the 
number of closed vortex-loops with a diameter given by 
the magnetic length in the problem. This number scales 
with L z at T z , and thus in the thermodynamic limit there 
are infinitely many such vortex-loops per flux-linc at the 
temperature where the FLL melts. 
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FIG. The distribution function D(p)/L 3 of thermally 
excited closed vortex loops as a function of perimeter p, 
for / = 1/32, r=3, L=48 and various temperatures. For 
T < Tbc2 — 1.8, D(p) is best fit to an exponential decay. For 
T > Tbc2, D(p) is best fit to an algebraic decay, indicating a 
vortex loop blow out at Tbc2 ^ T z > T m . 



D. Specific heat 



vant for discussing the low-field experiments of Zeldov et 
aln. Therefore, our simulations do not address the issue 
of the fate of the zero-field loop-transition in asymptoti- 
cally low magnetic fields, of the order of 100G and below. 
Neither do the simulations address how this finite-field 
counterpart of the zero-field vortex-loop transition inter- 
acts with the melting transition of the FLL and with 
the loss of phase-coherence along the field direction, in 
this Jpw-field regime relevant for discussing the results in 
Hef.fi 

We next turn to a discussion of the anomalies at the two 
lower temperatures T m and T z , and base our discussion 
on the lower panel of Fig. [?]. It shows the specific heat 
for the filling / = 1/32 and anisotropy V = 3, for various 
system sizes L = 16, 32, 48. The feature in the specific 
heat at T = T m is clearly associated with FLL melt- 
ing. One notable feature in the specific heat anomaly at 
T = T m is that its peak scales_as L 3 , characteristic of 
a first order melting transitions. This follows from the 
fact that a first order phase-transition is generally char- 
acterized by coexistence of two phases at the transition: 
One low-energy ordered phase and one high-energy disor- 
dered phase. Thus there is a discontinuity in the internal 
energy of the system at the transition temperature, and 
a delta-function peak in the specific heat. On a finite 
system the delta-function peak is converted to a peak of 
order L d , where L is the linear dimension of the system, 
and d is its dimensionality. In fact, the coefficient of the 
L d term is the discontinuity in the entropy of the system, 
at the transitions . so for the 3D case we obtain 



Our results for the specific heat are shown in the bot- 
tom panel of Fig|j| for the parameters / = 1/32, T = 1 
and L = 16,32,48,64,80,96. Similar results for the 
anisotropic case T = 3 are shown in the bottom panel 
of Fig. [?]. It is clear that that the lowest anomaly in 
the specific heat correlates with the melting of the FLL. 
There is also a broad feature at a much higher tempera- 
ture, Tbc2 ~ 1-9, which is roughly equal to the tempera- 
ture at which we see a sharp peak in the specific heat in 
the zero-field case. 

The feature at Tb C 2 is the remnant of the zero-field 
anomaly on the specific heat, previously shown in Fig|]. 
This broad peak in the specific heat is associated with 
the upper critical field H C 2 (T) , and our results show that 
in the extreme type-II Villain-model, H C 2 (T) line is very 
steep close to the zero- field T c . As discussed in the pre- 
vious section, Tb c i is also close to the temperature at 
where the distribution function D(p) of closed vortex- 
loops changes behavior for exponential to algebraic decay. 
Hence, we conclude that in the field-regime correspond- 
ing to / = 1/32, the loop-transition appears close to the 
mean- field 7J C 2(T)-line. 

We caution the reader that / = 1/32 is not a particularly 
low magnetic field. If we estimate it for YBCO using the 
method of Ref.El, it corresponds to a magnetic field pf 
the order of IT. As recently emphasized by TesanovicEiJ 
and Nguyen et aZ.til such magnetic fields may not be rele- 
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Thus we may use finite-size scaling of the specific heat to 
extract AS, or equivalently, the latent heat of the melting 
transition. We find 

AS ~ 0.00 /cs/vortex per layer, T = 1.0 

AS w 0.03 /c B /vortex per layer, T = 2.0 

AS w 0.05 fcs/vortex per layer, T = 2.5 

AS w 0.10 fc B /vortex per layer, T = 3.0 

The results for a finite-size scaling of the specific heat 
per site C and the entropy discontinuity at the melting 
transition, are shown in Fig. [Io|. AS* is seen to increase 
rapidly with T. This is expected on general grounds, 
since the flux-line liquid in a very anistropic supercon- 
ductor is expected to exhibit more disorder than in an 
isotropic case due to the more flexible nature of individ- 
ual flux lines in layered compounds. In our opinion, a 
reduction of AS with increasing T would be unphysical, 
at least in situations where the superconductor can be 
viewed as an anisotropic continuum. The 2D case is in 
some sense a singular limit, as discussed in Section II. D. 
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FIG. [l^. Top panel: The specific heat-maximum at the 
flux-line lattice melting transition, as a function of system 
size, for a given mass-anisotropy F = 2.5. Such plots may 
be used to extract the entropy discontinuity at the melting 
transition, as explained in the text. Bottom panel: Entropy 
discontinuity at the melting transition as a function of the 
mass-anisotropy. The first order character of the flux-line 
lattice melting is seen to increase rapidly with increasing V. 



Although it is expected on general grounds that the 
melting of the FLL is first order also in the isotropic 
case, it appears to be too small to detect in our simula- 
tions. However, by extrapolation we may easily extract 
entropy jumps of the correct order of magnitude seen 
in calorimetric data on YBCO with r w 83, since AS 
is expected to grow with further increase in T. If we 
extrapolate our results to T = 7, we obtain the value 
AS* = 0.35fcs/ per vortex per layer, in very good agree- 
ment with experimental data on YBCO. This estimate 
is again in surprisingly good agreement with, and only 
slightly larger than, the result obtainedjby Hetzel el. for 
the uniformly frustrated 3DXY-modeB with a consider- 



ably more sophisticated technique, but where the effect 
of anisotropy was not fully accounted for. 
There is also a weak specific heat anomaly at T = T z , 
associated with the loss of phase-coherence in the BCS 
order parameter along the field direction. What this 
specific heat anomaly conceivably could show, is that 
there is a phase-transition inside the flux-line liquid 
phase, from a low-temperature flux-line liquid phase, to 
a high-temperature ihrx-line liquid phase, as suggested 
by Feigel'man et aZ.E3. In RefEj, the low-temperature 
phase is suggested to correspond to a flux-line liquid with 
no entanglement of flux-lines, while the high-temperature 
flux- line liquid phase is suggested to correspond to a flux- 
line liquid with-fintanglement. In the language of the 2D 
boson-analogyo, the former would correspond to a nor- 
mal Bose-liquid, while the latter would correspond to a 
superfluid Bose-liquid, the two being separated by a gen- 
uine phase-transition. 

Note however, that the temperature T z goes down with 
increasing system size, apparently with no sign of satu- 
ration, approaching T m from above. What this strongly 
indicates, is that the portion of the phase-diagram with 
a flux-line liquid phase with an apparent intact phase- 
coherence in the BCS order parameter across the sample 
parallel to the magnetic field, will vanish in the thermo- 
dynamic limit. 

The question remains as to what the character of the flux- 
line liquid phase with no phase coherence along the mag- 
netic field, is. We have addressed the issue of whether 
the transition at T z — > T m results in a flux-line liquid 
with well-defined flux lines, by considering the amount 
of flux-line cutting and intersectioning of flux-lines with 
closed vortex loops, that takes place below and above the 
temperature T z — > T m . 



E. Flux-line cutting 



A flux-line liquid with large amounts of flux-line cut- 
ting events, is not likely to be able to sustain a heav- 
ily entangled vortex configuration with well defined flux 
lines. The amount of flux- line intersectioning, p Eq. (Q) 
is shown in Fig. |H} As we see, p increases sharply from 
zero at T — T z , and continues to increase monotonically 
as a function of temperature. Flux-line cutting is an effi- 
cient way of disentangling flux lines, and the large values 
of p suggest that in the flux-line liquid phase, above T z , 
the flux- line liquid is incapable of sustaining an entangled 
configuration. Hence, the loss of phase-coherence along 
the direction of the magnetic field essentially is due to 
intersectioning between flux-lines, and between flux lines 
and vortex-loops, with associated massive flux-line re- 
combinations. Recall that T = T z is also the temperature 
at which the number of closed vortex-loops of diameter 
equal to the magnetic length in the problem starts scal- 
ing with L z , and it is natural to associate the increase in 
p with this limited proliferation of vortex- loops in a finite 
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field. Under such circumstances, a world-line picture of 
2D bosons appears unlikely to be a particularly useful 
analogy to the flux-line liquid system. 
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FIG. [Tl| The number of vortex-intersection events per flux- 
line p, for / = 1/32, r = 3, and L = 48. 



V. SUMMARY 

As discussed in section II. C, there can only be one 
single 3D phase-transition in the superconductor in zero 
field, regardless of the anisotropy as long as this is finite, 

r < oo. 

We have shown that in zero field, the superconductor- 
normal metal transition is due to a vortex-loop transi- 
tion analogous to that first suggested to. occur in neutral 
superfluids, such as He 4 , by OnsageiO. The transition 
is characterized by a vanishing vortex-loop line tension, 
giving a loop-distribution function which decays as r~ a 
at the transition, where r is the radius of the loop. Be- 
low T c , the decay of the loop-distribution function clearly 
appears to be faster-than- algebraic. If the anisotropy is 
increased the vortex-loop unbinding temperature is re- 
duced, but the transition remains 31?. Close to the tran- 
sition the system is isotropized due to an upward renor- 
malization of the interplane coupling as a consequence 
of a diverging superconducting coherence length at the 
transition. 

An important point is that in finite magnetic fields, the 
situation is qualitatively different. In this case, the mag- 
netic length of the vortex system, i.e. the average dis- 
tance between the flux-lines, cuts off tke,-upward renor- 
malization of the inter-plane coupUngl3'tlL Note that 
the magnetic length only sets a new length scale for 
the phase-coherence as long as the flux-lines actually 
fluctuate. In the ground state, with no fluctuations in 



the Ahrikosov vortex lattice, the zero-field arguments 
applyB. 

Therefore, the perfect Abrikosov vortex lattice has no 
effect on the vortex-loop blowout transition, since the 
phase- coherence length is essentially infinite for this case. 
Conseguently, for the case A = oo, one cannot have a 
vortex-loop blowout transition far below the flux-line lat- 
tice melting temperature. 

As soon as the flux lines start to fluctuate appreciably, 
which only happens very close to the melting transition 
due to the first order character of the transition, the re- 
sulting much smaller cutoff on the rcnormalization of the 
inter-plane coupling facilitates a vortex loop blowout vir- 
tually at the same temperature as the flux lines start 
to fluctuate. We emphasize that these statements apply 
to the case of total suppression of gauge-fluctuations, or 
A = oo. 

Even if vortex-loops do not exist on all length scales in 
the problem at the flux-line lattice melting temperature, 
they will seriously affect the flux-line liquid phase pro- 
vided that they exist on length scales up to the mag- 
netic length of the problem. From our simulations, we 
have found that this is always the case in the thermo- 
dynamic limit.: The number of such closed non- field in- 
duced closed vortex-loops scales with the thickness of the 
sample, L z , whereas the number of field-induced flux- 
lines obviously does not. Therefore there is an infinite 
number of closed vortex-loops with a diameter equal to 
the magnetic length, per flux-line. Hence, in the liquid 
phase flux lines cannot be considered as well defined en- 
tities. They are well-defined entities in the lattice case. 
When a flux-line lattice melts, the molten phase is an 
incoherent vortex-liquid characterized by loss of phase- 
coherence along the direction of the magnetic field. This 
loss of phase-coherence is associated with a proliferation 
of closed vortex loops and massive flux-line cutting and 
recombination. Such a vortex-liquid phase is unlikely to 
sustain heavily entangled vortex-configurations. 
Note added: After this work was completed, we received 
works from Hu, Miyashita, and TachikiEU, as well as from 
Koshelevo, that support the conclusion that T z = T m 
when A =p-po, at least iniJie field regime considered, i.e 
/= 1/25E3 and /= 1/36E3. 
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APPENDIX A: HELICITY MODULUS IN TERMS 
OF PHASE- VARIABLES 

In this appendix we derive the expression for the he- 
licity modulus (Eq. ^]), for the uniformly frustrated 
anisotropic Villain model. 

The effective Hamiltonian for the uniformly frustrated 
anisotropic Villain model is 



Using the definition 

X = 9(f + e v ) - 9(f) - A v (r) + (e v ■ e^S, 

dx _ « ~ 
86 v ' 

we can write in the following form 



H v {6'(f)} = J J2 V ^ 6 '^ 



72 



k B T 
Jo 



r,u—x.y,z 
( 00 



In- 



E 



exp 



~2k B T 



(x — 2irmy 



(Al) 



+ 



Here, Jo, A v (f) and a u are define d in the text, e„ is the 
unit vector for the i>-axis. In Eq. Al, 9'(f) is the phase 
of the complex superconducting order parameter. {9'(f)} 
denotes the functional of 9'(f). We now apply boundary 
conditions such that the phase across the system in the 
/2-direction is twisted by an amount L5. If [1 = z, for 
example, we have the following phase twist 

6'(x,y,z = L)-9'(x,y,z = 0) = L5. (A2) 

Now we define a set of new phase variables 

6(f) = 9'(f)-(f-e fl )5. 

If the phase variables 9'(f) obey the twisted boundary 
conditions in Eq. A2, the new phase variables 9(f) obey 
the following periodic boundary conditions 



L 3 k B T 

Jo 
L 3 

Jf J_ 

L 3 k B T 



1 (T, WP)} [^^^M e 



k R T 



y [£ 



Z(<5=0) 
r K,(x)(e„-e M ) 2 ]e »bT 



E FE 



Z(<5=0) 



Z(<5=0) 



which is Eq. || written in a more explicit form. The sum- 
mations over the configurations {9(fi)} are restricted to 
configurations satisfying periodic boundary conditions. 



APPENDIX B: HELICITY MODULUS IN TERMS 
OF VORTEX SEGMENTS CORRELATIONS 

In this appendix we derive the helicity modulus for 
the anisotropic LSM T /i (fc), expressecLin terms of vortex 
density-density correlation functionsEd, /1 = (x, y, z). 
In the continuum limit, the effective Hamiltonian for the 
anisotropic LSM can be written as 



(x, y,z = L)- 9(x, y, z = 0) = 0. 



(A3) 



H{vJ} 



Ji 



d 3 7 



v(f)- M ■v(r) + 2k\ 2 f(f) ■ f(f) , (Bl) 



The Hamiltonian for the uniformly frustrated anisotropic 
Villain model in terms of the new phase variables is 

H v {6,6} = Jo V »W+ el) - 9(f) - A u (f) + (e v ■ e^)S}. 



The partition function in terms of the new phase vari- 
ables is 



Z(S) 



mm 



the energy per unit length Ji = $g/167r 3 A^. Here, v(f) 
is the super-fluid velocity, 

v(f) = V0(r) - 1(f), 

6(f) is the phase of the complex superconducting order 
parameter, &oA(f)/2ir is the vector potential, $0 is the 

flux quantum. In Eq. Bl, M is the anisotropic mass ten- 
sor, 



M= 



where we only sum over configurations {9(f)} that sat- 
isfy periodic boundary conditions. The total free energy 
is 

F(6) = -k B T \nZ(5). 

The helicity modulus T p is defined as the second deriva- 
tive of the free energy with respect to a phase twist across 
the sample in the /i-direction. Thus, 



1 
1 







r 2 



describing uniaxial c-anisotropy. The anisotropy param- 
eter F = A c /A a , A Q and A c are explained in the text. 
In Eq. Bl, f(f) is the local density of the magnetic flux 
quanta, 



B(f) 
$o 



-V x A(f). 



1 



d 2 F(S) 
dS 2 



Her e, B (f) is the local magnetic induction. In evaluating 
Eq. Bl, the integration must be cut off at the core of the 
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vortex segments, so that the energy stays finite. 
The partition function Z is computed averaging over in- 
dependently fluctuating v(r) and f(r), subject to the 
constraint (f(r)) = fe z for a constant uniform average 
magnetic induction Be Zl where (...) denotes a thermal 
average, and is the unit vector along the /2-axis. 
Using the Fourier transform 

iJ(jfe) = f d 3 r e*^(r), 



the Hamiltonian in Eq. Bl can be written as 

H{v(k)J(k)} = ^^2H'(k), 



H'(k) = v(k)- M -v(-k) + 2n\ 2 f(k) ■ f(-k) 



(B2) 



where V is the volume of the system. An applied twist in 
the phase of the superconducting order parameter, along 
the /i-direction, as considered in the previous Appendix, 
corresponds to a change in the super-fluid velocity 

v(k) — * v(k) + 8v{k)e^. 

In terms of velocities, therefore, the helicity modulus 
T M (fc) is defined as the second derivative of the free en- 
ergy F = —ksTlnZ with respect to a change in the su- 
perfluid velocity along the //-direction. Hence, 



d 2 F 



d[5v(k)]d[5v(-k)] 



V 



M -e. 



Sv(k)—Sv( 

Ji 



k B TV 



v(k)- M -e M ][e M - M -v(-k)} 



(3) 



The subscript denotes the unperturbed system with 
no applied phase twist, Sv(k) = 0. In Eq. ||, we have 
used (v(k))o = (v(—k))o = 0, since the Hamiltonian 
Eq. ~ 



contains only quadratic terms like [v v (k)v v {— k)] 
(y = x, y, z). Because of the symmetry in which V# and 
A enter v(r), T M (fc) can also be interpreted as the lin- 
ear response coefficient of the supercurrent induced by a 
perturbation in the vector potentialEU, 

= -T,(k)6A,(k). 

There are three interesting helicity moduli to be consid- 
ered: (1) T x (ke y ), which is the energy cost corresponding 
to a compressional perturbation of the flux-line system. 
(2) T x (ke z ), which is the energy cost corresponding to a 
tilting perturbation of the flux- line system. (3) T z (ke x ), 
which is the energy cost corresponding to a shearing of 
the flux-line system. 

To find an expression for T^(fce !/ ) (fi v) expressed in 



r 



terms of vortex segments density correlations, we need to 
write the parts of H containing ke v in a diagonal form. 
Defining the vortex segment density 



n(r) 



1 

2^ 



v x we, 



we can write the suoer-fluid velocity in the following 
gauge-invariant formE3 



v(k) = 2m kx{k) 



k x [n(k) - f(k)} 
k 2 



(4) 



Here, x(r) is a smooth scalar function which describes the 
longitudinal part of v(k). The transverse part of v(k) is 
determined by V x v(r) = 27r[n(r) — /(r)]. 
Substituting Eq. || into H'(k) (Eq.[B^), we get the follow- 
ing diagonal form for H'(k) for the case of k = ke v 



H'(ke x ) 
4tt 2 

H'(ke y ) 
4tt 2 

H'{ke z 
4tt 2 



k 2 x(ke x )x(~ke x ) 



k 2 x(ke y )x(~key 



Ai{k)n y (ke x )n y {-ke x ) + Ai(k)n z (ke x )n z {~ke x ) + 
B x (k)^ y {ke x )^ y {~ke x ) + B 2 (k)Sf z (ke x )5f z (~ke x ) 



Ai(k)n x (ke y )n x (-ke y ) H 
Bi(k)Sf x (ke y )Sf x (-ke y ) 



A 2 (k)n z (ke y )n z (-ke y ) + 
+ B 2 (k)Sf z (ke y )Sf z (-ke y ), 



= (k/T) 2 x{ke z )x{-ke z ) + A 2 (k)n x (ke z )n x (-ke z ) + A 2 (k)n y (ke z )n y (-ke z ) + 

B 2 (k)Sf x (ke z )6f x (-ke z ) + B 2 (k)Sf y (ke z )Sf y (-ke z ). 



(5) 



Here. 



Ax(k) = 
Bi(k) = 



A 2 



l + T 2 A 2 g 2 ' 
l + r 2 A 2 g 2 



rv 



A 2 (k) = 
B 2 (k) = 



A 2 



1 + AV 
1 + X 2 q 2 
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In Eq. [|, 6f(k) — f{k) — f°(k) is the fluctuation of the magnetic flux density away from the value f°(k) minimizing 
the Hamiltonian for a given vortex segments configuration^. 



f°(k) 



h{k) 



(r 2 -l)[n(fc).q] A 2 <f 



i + \ 2 k 2 i + A 2 fc 2 + (r 2 -i)AV 



where q = k x z. To compute the partition function Z, we should sum over: (1) all smooth functions (2) all 

h(k) that satisfy k ■ n(k) = 0, and (3) all 8f(k) that satisfy k ■ 8f(k) = 0. The constraints (2) and (3) come from the 
restriction of no divergence in the vortex segments density and no divergence in the local magnetic field induction. 
Substituting Eq. |] into Eq. El, and using the the Hamiltonian Eq. ^ to evaluate the average over x(k) and Sf(k). For 
the case of k = ke v we obtain 



T^(ke u ) 



Ji 



X 2 k 2 



1 



4-k.Ji A 2 (n< T (fce„)n CT (-fce„))o 



v i + [i + v,(r 2 - \)]\ 2 k 2 V k B TV i + [i + <5 M , 2 (r 2 - i)]A 2 P 

(fx, v, a) are cyclic permutation of (x, y, z). The generalization of Eq. |to a lattice superconductor is 
Jo (A/d) 2 Q 2 / 47tJ (\/d) 2 (n CT (fce l/ )n CT (-fce I/ )) 



+ <W(r 2 -i)](A/d) 2 Q 2 



fc B TL 3 i + [n-^,,(r 2 -i)](A/d) 



2/Q2 



(6) 



(7) 



where Q M = 2sin(k fl d/2), Q 2 = ^2 Q 2 , fc M is the /i-component of k, and d is the lattice constant of the (simple cubic) 
underlying numerical lattice. 
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